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Abstract: Determination of site quality is a basic tool for proper selec¬ 
tion of locations and species, in management of forest plantations. 
Throughout the Caribbean studies of site quality are few and are ham¬ 
pered by statistical limitations, inappropriate growth models, and limited 
data. We fitted growth curves for dominant height to evaluate and clas¬ 
sify site quality of teak (Tectona grandis) plantations by using data from 
44 permanent sample plots established since 1990 in 34>2 years old teak 
plantations in the Colombian Caribbean region. We used Korf’s and von 
Bertalanffy’s models to fit curves as non-linear effects models. Both 
models, with a single random parameter, were considered as adequate for 
dominant height growth modelling, but Korf’s model was superior. The 
resulting curves were anamorphic and closely reflected high variability in 
site quality. Five site classes were clarified: at a base age of 12 years old, 
teak reached a mean dominant height of 24.8 m on the best sites, 9.8 m in 
the worst sites, and in the averages sites, 15.8-18.8 m. Using this model, 
we identified the best and the worst sites for teak plantations in the Car¬ 
ibbean region. This model proved a useful tool, not only for site quality 
evaluation, but also for improved teak plantation planning and manage¬ 
ment. 

Keywords: Tectona grandis; site index; growth modelling; permanent 
sample plots; Colombia 

Introduction 

Forest site quality refers to the sum of factors that affect the 
productive capacity of forests. Factors can be classified into three 
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groups: climatic, biotic, and edaphic (Carmean 1979). Site quality 
is defined as the potential production of timber in a particular site 
for a particular species or a specified forest type (Clutter et al. 
1983). The most commonly used method is direct sampling of the 
mean height of dominant trees in a stand, because in even-aged 
stands it is little affected by stand density or thinning (Clutter et al. 
1983). Site quality correlates closely with timber production in 
terms of volume or biomass (Clutter et al. 1983). According to this 
method, site index ( SI) is defined as a measure of site quality, 
based on the mean height of dominant and co-dominant trees of 
arbitrarily sampled age classes (Carmean 1979). The curves of SI 
that result from classifying and categorizing height growth curves 
can be fitted in several ways, but always through height-age 
coordinate pairs. 

Previous studies used the SI method for determining site qual¬ 
ity in teak (Tectona grandis L.f.) plantations, but most used data 
from temporary plots (Nunifu and Murchinson 1999; Henao 
1982; Keogh 1981 and 1982; Bermejo et al. 2004). Quite often, 
when Permanent Sampling Plots (PSPs) were used, the possible 
existence of polymorphism was not corroborated (Vaides et al. 
2004; Upadhyay et al. 2005; Jerez-Rico et al. 2011). Most re¬ 
searchers used the guide curve method (or proportional curves 
method), which generates anamorphic curves and uses the 
non-versatile Schumacher’s model (Nunifu and Murchinson 
1999; Henao 1982; Keogh 1982; Bermejo et al. 2004; Vaides et 
al. 2004; Jerez-Rico et al. 2011). Although this method has re¬ 
ceived little attention, it is important because the presence of 
polymorphic patterns can indicate two things. First, the height 
growth of a species is sensitive to silvicultural treatments. Sec¬ 
ond, silvicultural treatments are often continued after trees have 
reached heights where the treatments no longer have an impact. 
Therefore, one species may show both patterns of growth. 
Upadhyay et al. (2005) used both the guide curve and the differ¬ 
ence equation method with PSPs to develop site index curves in 
teak plantations in India using the Hossfeld IV growth model. 
They found the polymorphic difference equation method to be 
superior. From research on teak management regimes, it can be 
inferred that growth in height is independent of stand density. A 
15 year stand of teak in Nigeria received three treatments: un¬ 
thinned (2,200 trees- ha' 1 ), thinned to 760 trees- ha' 1 , and 
thinned to 395 trees- ha' 1 (Lowe 1976). When the stand was 20 
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years old, the fastest tree growth was unrelated to the treatments. 
Vincent et al. (2000) carried out a thinning experiment in Barinas, 
Venezuela. In a 13-year-old teak plantation, they established 
several spacing regimes from 2 m x 2 m to 4 m x 4 m. At 18.7 
years, there was no relationship between dominant tree height 
and thinning regime. Jerez-Rico et al. (2011), working in the 
western plains of Venezuela, examined data from permanent and 
temporary plots of teak that included more than 30 years of 
measurements. Initial spacing varied from 2mx 2 m to 4 m x 4 
m, and density at sampling ranged from 200 to 2,400 trees- ha' 1 . 
They found that dominant height was little affected by tree den¬ 
sity. On the Peninsula of Nicoya, Costa Rica, Chaves et al. 
(2003) established a thinning experiment in a seven year old teak 
plantation with the following treatments: 15 m 2 - ha' 1 , 17 
m 2 - ha' 1 , 19 m 2 - ha" 1 , 21 m 2 - ha" 1 , and 25 m 2 - ha" 1 . The experi¬ 
ment continued until trees reached 20 years in age. Growth in the 
dominant tree height was not affected by treatments. Jerez-Rico 
et al. (2011) used mixed models to study the growth of teak. 

The first approaches in Colombia for studying the growth of 
teak were from Echeverri (1968) and Rodriguez (1968), but they 
did not attempt to develop site index curves. Henao (1982), in a 
plantation of the Department of Cordoba and on the basis of 
average tree height from temporary sample plots, found no sig¬ 
nificant quality differences between sites. A recent and continu¬ 
ing study of SI is the site classification chart for the Caribbean, 
Central America, Venezuela, and Colombia (Keogh 1982). In our 
study, the data for Colombia were taken from the Venegas report 
in 1977 “Reply to teak questionnaire”, performed by the FAO 
COL/74/005 Project. In recent decades, various timber compa¬ 
nies planted teak not only in Colombia but also in other Carib¬ 
bean countries. However, the foresters from the Caribbean use 
few teak site index studies to properly manage this species. In 
Colombia, we used the result of Keogh (1982) that was based on 
limited data and used outdated methods. Our study aims to use 
modem statistical techniques and widely accepted growth mod¬ 
els to develop a family of site index curves for teak from the 


Colombian Caribbean. Our results are applicable to all Caribbean 
teak plantations. 

Materials and methods 

Study area 

In Colombia, and particularly in the Cordoba Department, Co¬ 
lombian Caribbean region, teak has been planted for over 70 
years (Rodriguez 1968; Echeverri 1968; Keogh 1981). According 
to the Holdridge (1982) classification scheme, the area falls un¬ 
der the lowlands monsoonal association of the humid tropical 
forest life zone. The rainfall regime is unimodal and at the most 
representative weather station (IGAC 1978) averages 2,479 
mm/yr. Approximately, 30% of rain falls between August and 
September and less than 60 mm per month falls during January 
and Febmary. The annual mean temperature is 27°C, with daily 
amplitude of 10°C. Relative humidity varies from 75% to 84%. 
All teak plantations are on colluvial-alluvial soils with flu¬ 
vial-lacustrine deposits used for livestock grazing before the 
plantations were established (IGAC 1989). 

Data source 

In 1990, 20 Permanent Sampling Plots (PSPs) were established in 
the teak plantations of the Cordoba Department, and the number 
of PSPs was increased in subsequent years (Fig. 1). The age of the 
PSPs ranged from 3 to 22 years. In 2011, there were 44 PSPs. The 
PSPs had two plot sizes. The plot size of 23 PSPs was 600 m 2 (i.e.: 
30 m x 20 m) and the plot size of the remaining 21 PSPs was 
1,000 m 2 (i.e.: 40 m x 25 m). Diameters at breast height (1.3 m 
above ground) and dominant height (mean tree height of the 100 
tallest trees per hectare) were measured annually. Dominant 
height was used as an indicator variable for site quality. 



Fig. 1 Study area. The dots show the location of the permanent plots. 

All the PSPs were subjected to the same management regime. spacing of 2.5 m x 2.5 m), manual weeding was carried out three 

The plantation stocking density was 1,600 seedlings- ha" 1 (i.e., a times during the first two years and once each year in subsequent 
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years until the end of the rotation. Pruning was carried out in 
years five and nine. Basal area was maintained close to 26 
m 2 - ha" 1 by two thinnings between years 7-9 and 1243. 

Modelling 

In order to fit mean dominant height growth functions, two mod¬ 
els used in silviculture were evaluated: von Bertalanffy (or 
Chapman-Richards’s, Eq. 1), and Korf (Eq. 2) (Kiviste et al. 
2002 ). 

Hd = A [l - exp (- A Of 2 (i) 


Hd = A exp 



where, Hd is the dominant height (m), A is the asymptote (m) or 
maximum value reached by Hd , fi\ and p 2 are the unknown pa¬ 
rameters, t is the age (year) corresponding to each//A, and exp() is 
the exponential operator (Euler’s constant). 

As all data were collected in PSPs, this is a typical case of re¬ 
peated measures over time (or longitudinal data) without inde¬ 
pendence among intra-plot measurements. Therefore, there can be 
serious problems in autocorrelation of errors, when the parameters 
are estimated by conventional least squares methods (linear or 
non-linear). Besides the problem of autocorrelation, growth and 
yield data from permanent plots usually exhibit heteroscedasticity 
(Gregoire 1987). 

A non-linear mixed-effect model was used to derive reliable 
estimators of the growth model parameters. It enabled modelling 
of the intra-individual covariance structure, assuming that (usu¬ 
ally one or two) individual-independent, small-dimension latent 
random-effects vectors existed in the model (SAS 1999). To apply 
this method, a re-parameterization of models (Eq. 1 and Eq. 2) 
was made (Fang and Bailey 2001), by changing the value of the 
asymptote A by an unknown parameter fi\. The value of fi\ cor¬ 
responds to the expected value of Hd when t = t 0 . Then Eq. 1 and 
Eq. 2 can be represented in the following way: 


Hd - A 


1-exp {-P 2 t) 



l-exp(-/? 2 0 



Hd = A exp 



1 






Where, [i\ = (p + b, unknown parameter corresponding to the SI 
value for a base age of t 0 = 12 years. A mixed parameter is con¬ 
sidered fixed for both models (i.e., it has a fixed (<p) and a random 
(b) part), and p 2 and p 2 are unknown parameters. 

Using this model, polymorphic and anamorphic growth curves 


can be obtained. Several authors (Cieszewski and Bella 1989; 
Cieszewski and Bailey 2000; Fang and Bailey 2001) examined 
the properties of these models. Eq. 3 and Eq. 4 are fitted follow¬ 
ing a single random effect in the Pi parameter, and each model is 
evaluated considering three variance structures: constant vari¬ 
ance, variance as a potential function of the mean, and variance 
as an exponential function of the mean, for a total of six models. 
Each model is evaluated in terms of the statistical significance of 
the parameters, and is compared through both the Akaike infor¬ 
mation criterion (AIC) and the Bayes information criterion (BIC). 
The selected model undergoes a diagnostic of residuals, as pro¬ 
posed by Fang and Bailey (2001). This modelling phase was 
performed using the NLM Procedure in SAS System for Win¬ 
dows V. 8 (SAS 1999). 

Site classes 

After selecting the best model, we estimated the SI for the average 
site (random parameter = 0). We then identified sites above or 
below the average by considering the variation of the random 
parameter. We then ranked sites according to the intervals of four 
standard deviations (± 2 standard deviations around the average 
site). 

To place each plot in its corresponding site class, the SI for 
each plot measurement was estimated (by isolating the Pi pa¬ 
rameter of the selected model) and the value was placed in the 
corresponding category. This procedure is also useful to evaluate 
the goodness of fit, because all of the measurements from the 
same plot should be placed in the same site class. Another way in 
which site classes can be assigned is by plotting on a coordinate 
diagram the curves that outline each site class, and then plotting 
the estimations of Hd for each plot. 


Results 

Exploratory analysis of the database 

Scatter plots of the dominant height for some of the PSPs are 
depicted in Fig. 2. Not all 44 PSPs are shown in Fig. 2, because 
many PSPs overlapped, and it was difficult to correctly interpret 
the scattering. However, the scatter plots show high variability in 
the growth of teak in the study area, reaching dominant heights at 
16 years ranging from 11-27 m. It is important to highlight the 
heteroscedastic nature of the scattered plots in Fig. 2 because 
variability increased with age. In most but not all cases, domi¬ 
nant trees sampled in successive sampling periods were the same. 
Therefore, there was temporal autocorrelation in the series of 
dominant height, another problem for the regression analysis. 

Modelling 

The parameters of Eq. 3 and 4, evaluated on each one of the 
variance structures, were statistically significant with a confi¬ 
dence level of 95% (Table 1). The values of the two selection 
criteria for each model are shown in Table 1 under the three 
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variance structures. 



Fig. 2 Mean dominant height versus age for 14 out of 20 PSPs 
measured nine or more times. Points with equal shape belong to the 
same plot 


Though all of the models were statistically acceptable, the 
models with variance as an exponential function of the mean were 
better (3.3 and 4.3). Of two values, Korf s model (4.3) was the 
most appropriate, because its Bayesian’s criterion was slightly 
lower than for von Bertalanffy’s model (Table 1) with the vari¬ 
ance of Hd as an exponential function of its mean, and a base age 
of 12 years (Eq. 5). There were no evident autocorrelation (Fig. 3), 
heteroscedasticity, or lack of fitness problems in relation to ob¬ 
servation order (observation number or the predicted values). The 
lack of outliers in Fig. 3 suggests a normal error distribution. The 
statistical results show that the non-linear mixed-effect model 
corrected the residual problems inherent in longitudinal data. It 
should be noted that von Bertalanffy’s model also shows similar 
characteristics of residuals (Fig. 3) 


Hd = 17.28 exp 
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Fig. 3 Residuals of the dominant height (Hd) as a function of both, observation order and estimated value (Hd est.) 


Table 1. Results of the information criteria: Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) a 




3 (von Bertalanffy) 



4 (Korf) 


Model 

3.1 

3.2 

3.3 

4.1 

4.2 

4.3 

Variance structure 

Constant 

Potential 

Exponential 

Constant 

Potential 

Exponential 

* 

AIC 

586.5 

583.4 

579.1 

578.5 

575.2 

570.7 

* 

BIC 

591.2 

588.1 

583.8 

583.3 

579.9 

475.4 


a The lower the value, the better the model, von Bertalanffy’s and Korf s growth models were similar and statistically acceptable. Korf s model (4.3) was 
the most appropriate because it achieved the minimum Bayesian criteria. 

Eq. 4 can be also written with site SI as a dependent variable of years, and by replacing the values of the other parameters pre¬ 
age (t) and Hd. This is accomplished by isolating the parameter sented on Eq. 5. In this way, Eq. 6 can be obtained, 

which, as mentioned above, corresponds to SI at base age of 12 
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Site index classes 

When intervals of 12 standard deviations are considered around 
the average site, five intervals are needed to cover the complete 
range of Hd. SI limits for each category are shown in Table 2 and 
the corresponding SI curves are shown in Fig. 4. There was great 
variability in the site indexes. For trees twelve years old, the site 
index curves covered plantations from approximately 9.8 m up to 
24.7 m of dominant height. 


Table 2. Limiting values for site index classes for teak in Colombia 
separated by intervals of ±2 standard deviations around the average 
site 


Site index 

classes 

SI 

' Number of PSPs 

Lower limit (m) 

Upper limit (m) 

I 

21.7677 

24.7617 

4 

II 

18.7737 

21.7677 

12 

III 

15.7797 

18.7737 

17 

IV 

12.7857 

15.7797 

8 

V 

9.7917 

12.7857 

3 
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Fig. 4 Site index curves for teak (12 years base age). Dotted lines are 
the SI curves and continuous lines are limits of each SI classes (I to V) 


Discussion 

Estimated values of SI were similar for different ages of the same 
PSP. Fig. 2 and Fig. 4 show high variability of Sis. Dominant 
height followed an anamorphic curve, and this was confirmed by 
the satisfactory fit of this kind of model. As discussed in the 
introduction, there is abundant evidence that height increase in 
teak is poorly related with stand density. This result differed from 
the findings of Fang and Bailey (2001), who researched slash pine, 
and considered all parameters as mixed due to the high variability 
of silvicultural treatments. They found height increments that fit 


polymorphic curves. In our study, only a random parameter was 
necessary due to the anamorphic nature of the curves. 

Height growth of the dominant teak trees was well represented 
by Korf s model as modified (Fang and Bailey 2001). Moreover, 
this model eliminated the incompatibility problem between tree 
height growth and SI described by Curtis et al. (1973). Another 
advantage of the model is that it does not vary with changes in 
base age (Bailey and Clutter 1974). This means that the SI pre¬ 
dicted by the model is independent of the reference age (Fang and 
Bailey 2001). 

In our plantations, 39% of the PSPs were categorized as the 
average site class III. The proportion of PSPs with SI classes 
above the average (Classes I and II) was 35%, while 26% of PSPs 
were below the average site class (Classes IV and V). The shape 
of SI curves resembled those of Henao (1982) (Fig. 5a). In the 
value range, they resemble more those of Keogh (1982) (Fig. 5b). 
This last result is interesting, because in the Colombian Caribbean 
all of the SI classes reported by Keogh (1982) for Central America, 
the Caribbean, Venezuela, and Colombia were represented. The 
site ranges, presented by Miller (1969) for Trinidad, by Bermejo 
et al. (2004) and Perez and Kanninen (2005) for Costa Rica, by 
Jerez-Rico et al. (2011) for Venezuela and by Vaides et al. (2004) 
for Guatemala, also fell within the SI curves developed in this 
study. However, top height in Ghana at 12 years varied between 8 
m on the poorest sites and 15 m on the best sites (Nunifu and 
Murchinson 1999), figures substantially lower than in the Carib¬ 
bean. 

Both Henao (1982) and Keogh (1982) in Colombia used 
Schumacher’s model (Schumacher 1939), which is a simplifica¬ 
tion of Korf s model (Eq. 2) because the exponent of t is equal 
to one. This simplification reduces the versatility of Korf s model, 
compelling the function to reach a change in concavity (Point of 
inflexion in t = pi/2, and Hd = 0.135^4) at exactly 13.5% of the 
asymptotic value without relation with the recognized explosive 
initial height growth of teak. This is why the work of Henao and 
Keogh underestimated the initial growth and the SI of teak (Fig. 
5a and 5b). 

Upadhyay et al. (2005) used the Hosfeldt IV model (Kiviste et 
al. 2002) for top height SI curves at 25 years base age for teak in 
India according to 150 PSPs established in all teak plantations in 
the country, representing the top height from 4 to 93 years. Similar 
sites were estimated by these curves at their base age, an approach 
different to that used in this study. At 25 years the poorest sites 
were the same as in our study {Hd = 11 m) and the best sites were 
similar (29 m in our study versus 28 m in the Indian study). 
However, the slope of SI curves in our study began higher and 
then declined, while the SI curves of the Indian study maintain a 
more stable slope. Both mean SI curves intersected at t —32 
years. Therefore, after this year the Indian SI curves predicted 
higher top height than dominant height as found in our study. The 
opposite occurred before 32 years of age (Fig. 6). In spite of the 
large database in the Indian study, they did not measure, but rather 
estimated the top height of the trees by using an allometric rela¬ 
tionship between the quadratic mean diameter of ten trees with 
larger diameter per plot versus the average height of these trees 
(Upadhyay et al. 2005). Because trees with the largest diameter 
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are not necessarily the tallest, this procedure tends to overestimate 
the top height of the older trees with larger diameters. This is due 
to the fact that old trees continue growing in diameter but do not 
necessarily grow more in height. The same pattern results when SI 
curves of Upadhyay et al. (2005) are compared with many of the 
published teak SI curves (Miller 1969; Henao 1982; Keogh 1982; 
Malende and Temu 1990; Nunifu and Murchison 1999; Bermejo 
et al. 2004; Mora and Meza 2004; Perez and Kanninen 2005; 
Jerez-Rico et al. 2011). 




Fig. 5 Comparison of SI curves in the study (solid lines) with (a) 
Henao (1982) and (b) Keogh (1980) SI curves (dotted lines) 



Fig. 6 Comparison between the mean SI curve in t the study (thick 
dark line) and Upadhyay et al. (2005) mean SI curve (thin clear line) 

Growth of dominant tree height for teak in Colombia reflects 
site quality. The array of sites yields an anamorphic system of 
curves, allowing site classification for this species. In the Co¬ 


lombian Caribbean region, most previously reported ranges of SI 
for teak throughout the tropics are represented. 

Potential users of this study for the site index classes of a teak 
plantation must use the same criterion for dominant height as 
used in this research for modelling SI curves. In this way, with 
the knowledge of the dominant height (i.e. the mean height of 
100 tallest trees per hectare ( Hd) or its estimation through per¬ 
manent or temporary sample plots) and the corresponding age, 
the coordinate pair ( Hd , age) can be plotted as in Fig. 4. Thus, 
the site class to which the measures correspond can be identified. 
In order to calculate the SI, that is, the exact value of Hd ex¬ 
pected at 12 years, Eq. 6 must be used. This value can be catego¬ 
rized afterwards, when confronted with the limit values of each 
SI class presented in Table 2. 
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